
import MatrixIO
import cern
import jbcl,nix
from matrix import Matrix,size
import java,javax,jarray

from affymetrix.fusion.cdf import *
from affymetrix.fusion.cel import *

## this is a straight java port of RMAexpress c/c++ code. Any errors are mine.
## [sgalbrai@cs.ucla.edu]

def qnorm(X):
	# maybe I can write to a file, then load, then do the qnormalization.
	qn = nix.affy.QuantileNormalization(0, False, True)
	t = X.toArray()
	qn.normalize(t)
	X.assign(t)

def get_mu(X):
	PM = cern.colt.list.DoubleArrayList(X)
	stat = cern.jet.stat.Descriptive
	return stat.mean(PM)

def get_sd(X):
	PM = cern.colt.list.DoubleArrayList(X)
	stat = cern.jet.stat.Descriptive
	ss = stat.sumOfSquares(PM)
	sum = stat.sum(PM)
	var = stat.variance(PM.size(),sum,ss)
	sd = stat.standardDeviation(var)
	return sd

def get_alpha(PM,PMmax,rows):  # PM is a column vector
	alpha = 0.0;
	tmpsum = 0.0;
	numtop = 0;
	for i in range(0,rows):
		if (PM[i] > PMmax):
			tmpsum = tmpsum + (PM[i] - PMmax);
			numtop = numtop + 1;
	if tmpsum==0.0:
		alpha=0
	else:
		alpha = numtop/tmpsum;
	return alpha ;


def kernelDensity(PM):
	kernel = jbcl.calc.statistics.kernels.NormalKernel()
	#bd = kernel.findOptimalBandwidth(PM)
	bd = 1.0
	kd = jbcl.calc.statistics.OneDimKDE(PM, bd, kernel)
	return kd

def get_param(v,column,kd):
	d=0
	idx=0
	for x in xrange(0,len(v)):
		mt = kd.evaluate(x)
		if d<mt:
			idx=x;
			d=mt
		else:
			pass

	PMmax = idx;
	sd = get_sd(v);
	alpha = get_alpha(v,PMmax,len(v));
	param = [alpha,PMmax,sd]
	return param

def Phi(x):
	seed = java.util.Date()
	return cern.jet.random.Normal(0,1,cern.jet.random.engine.MersenneTwister(seed)).cdf(x)

def phi(x):
	seed = java.util.Date()
	return cern.jet.random.Normal(0,1,cern.jet.random.engine.MersenneTwister(seed)).pdf(x)
 
def bkgCorrect(M):
	rows=M.rows()
	columns=M.columns()
	kd = kernelDensity(M.like1D(rows*columns).toArray())
	##kd = kernelDensity(M.elements)  #need to fix this as it is slowing things down considerably.
#	print rows
	for j in xrange(0,rows):
		v = M.viewRow(j).toArray()
		###kd = kernelDensity(v)  #need to fix this as it is slowing things down considerably.
		param = get_param(v,j,kd)
		### IF I CAN DO THIS BY ROW IT WILL BE FASTER!!
		for i in xrange(0,columns):
			a = M.getQuick(j,i) - param[1] - param[0]*param[2]*param[2];
			r = a + param[2] * phi(a/param[2])/Phi(a/param[2]);
			M.setQuick(j,i,r)
		print "processed experiment", j
	kd=None
	return M
	

def median_polish(data,rows,cols,nprobes):

	# get each row median
	# get each col median
	# subtract from each row/column
	pass

 
  

class ExpressionSet:
	cdf = None; 
	CDF_FILE = 'HG-U133A.cdf'
	CEL_FILES = []

	def __init__(self,directory):
		dir = java.io.File(directory)
		self.CEL_FILES = dir.listFiles()

	def processCEL(self):
		# read first cell file, get dimensions
		###debug self.CEL_FILES = self.CEL_FILES[0:3]
		EXPS = len(self.CEL_FILES)
		cel = FusionCELData()
		M=None
		j=0
		for f in self.CEL_FILES:
			cel.setFileName(f.getAbsoluteFile().toString())
			cel.read()
			ROWS = cel.rows
			COLS = cel.cols	
			if M is None:
				# do the matrix allocation here.
				M = cern.colt.matrix.impl.DenseDoubleMatrix2D(EXPS,ROWS*COLS)
			for x in xrange(0,ROWS):
				for y in xrange(0,COLS):
					idx=cel.xyToIndex(x,y)			
					M.setQuick(j,idx,cel.getIntensity(x,y))
			j=j+1

		M = bkgCorrect(M)
		qnorm(M)
		F = cern.colt.matrix.linalg.Algebra()
		print F.transpose(M)
		
	def readcdf(self):
		self.cdf = FusionCDFData();
		self.cdf.setFileName(self.CDF_FILE);
		self.cdf.read()
	
	def printCDF(self):
		for i in range(0,S):
			name = self.cdf.getProbeSetName(i)
			print name
			self.cdf.getProbeSetInformation(i, set);
			assert( set.getProbeSetType(), FusionGeneChipProbeSetType.ExpressionProbeSetType );
			GROUPS = set.getNumGroups();
			for g in range(0,GROUPS):
				set.getGroup(g, group);
				print "\t", group.getName()
				for y in range(0,set.getNumCells()):
					group.getCell(y, probe);
					print "\t\t", probe.getX(), probe.getY();

	def RMA(self):
		Q = bkgCorrect(self.PM)
		S = qnorm(Q)
		E = summarize(exprSet)


class AffyApplet(javax.swing.JApplet):
# provide the GUI and interaction
	def __init__(self):
		pass

	def loadCDF(self):
		pass

	def getLocalDirectory(self):
		pass



### test code ###

#filename = 'c:\\dist\\Etest.txt'
#[X,Xl,Cl] = MatrixIO.readmatrix(filename,None,False)
#w = X[:,5]._M.viewColumn(0).toArray()
#Q = bkgCorrect(X)

x = ExpressionSet('cel/')
x.processCEL()


#header = x.cdf.getHeader();
#S = header.getNumProbeSets();
#set= FusionCDFProbeSetInformation()
#group = FusionCDFProbeGroupInformation()
#probe = FusionCDFProbeInformation()
## I only need to worry about the CDF file at the end.. to do the summarization.
## just go through the CEL files




### read each cel file into a matrix as  [(x,y), experiment]
